Data visualizations

Point estimates

results_scaled_group %>% 
  filter(adj == "Age, gender and weight change") %>% 
  select(variable, estimate, conf.low, conf.high, p.value, p.value.group, fdr.group) %>% 
  left_join(annotation$nightingale, by = c("variable" = "name.short")) %>% 
  filter(type == "Particle concentration") %>% 
  
  # Initiate ggplot call, and put "variable" and "estimate" on the x and y axes, respectively
  ggplot(data = ., mapping = aes(x = variable, y = estimate)) + 
  
  # Add a geom_point layer (a geometric point object)
  geom_point()

This is a start. We have our variables and our estimates. But what we really want is to also visualize the uncertainty around those estimates. Which geometrical object should we therefore switch to?

Point estimates with uncertainty

results_scaled_group %>% 
  filter(adj == "Age, gender and weight change") %>% 
  select(variable, estimate, conf.low, conf.high, p.value, p.value.group, fdr.group) %>% 
  left_join(annotation$nightingale, by = c("variable" = "name.short")) %>% 
  filter(type == "Particle concentration") %>% 
  
  # Before this point: dplyr specifics (data wrangling)
  
  # Below this point: ggplot2 specifics (plotting) 
  # Note also that I use "+" and not "%>%" after this point
  
  ggplot(mapping = aes(x = variable, y = estimate)) + 
  
  # Switch to geom_pointrange
  geom_pointrange()

But there is something wrong! What does the error message say?

Ah, so we might need to add the lower and upper CI boundaries. What were thir names again?

names(results_scaled_group)

##  [1] "adj"           "variable"      "term"          "estimate"     
##  [5] "std.error"     "statistic"     "p.value"       "conf.low"     
##  [9] "conf.high"     "r.squared"     "adj.r.squared" "sigma"        
## [13] "statistic1"    "p.value1"      "df"            "logLik"       
## [17] "AIC"           "BIC"           "deviance"      "df.residual"  
## [21] "fdr"           "p.value.group" "fdr.group"

They are named “conf.low” and “conf.high”! Perfect. Let’s add those to the geom_pointrange aesthetics.

results_scaled_group %>% 
  filter(adj == "Age, gender and weight change") %>% 
  select(variable, estimate, conf.low, conf.high, p.value, p.value.group, fdr.group) %>% 
  left_join(annotation$nightingale, by = c("variable" = "name.short")) %>% 
  filter(type == "Particle concentration") %>% 
  
  # Before this point: dplyr specifics (data wrangling)
  
  # Below this point: ggplot2 specifics 
  # Note also that I use "+" and not "%>%" after this point
  
  ggplot(mapping = aes(x = variable, y = estimate)) + 
  
  # Add a geom_pointrange layer with lower and upper CIs
  geom_pointrange(mapping = aes(ymin = conf.low, ymax = conf.high))

But what we really want is to display the variables on the y axis. Flip the axes using “coord_flip”.

Flip coordinates

results_scaled_group %>% 
  filter(adj == "Age, gender and weight change") %>% 
  select(variable, estimate, conf.low, conf.high, p.value, p.value.group, fdr.group) %>% 
  left_join(annotation$nightingale, by = c("variable" = "name.short")) %>% 
  filter(type == "Particle concentration") %>% 
  
  # Before this point: dplyr specifics (data wrangling)
  
  # Below this point: ggplot2 specifics 
  # Note also that I use "+" and not "%>%" after this point
  
  ggplot(aes(x = variable, y = estimate)) + 
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) + 
  
  # Flip the coordinates
  coord_flip()

Nice. But now it becomes apparent that the variable names are ordered alphabetically. That’s really not helpful at all. Use the fct_reorder function to reorder the variables by the “name.order” column.

Reorder variables

results_scaled_group %>% 
  filter(adj == "Age, gender and weight change") %>% 
  select(variable, estimate, conf.low, conf.high, p.value, p.value.group, fdr.group) %>% 
  left_join(annotation$nightingale, by = c("variable" = "name.short")) %>% 
  filter(type == "Particle concentration") %>% 
  
  # Reorder the variables 
  mutate(variable = fct_reorder(variable, name.order)) %>% 
  
  # Before this point: dplyr specifics (data wrangling)
  
  # Below this point: ggplot2 specifics 
  # Note also that I use "+" and not "%>%" after this point
  
  ggplot(aes(x = variable, y = estimate)) + 
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) + 
  coord_flip()

Better - but it feels natural to keep the largest size high up on the plot, and the smaller sizes lower on the plot. So simply add a minus sign to the ordering to switch it around.

results_scaled_group %>% 
  filter(adj == "Age, gender and weight change") %>% 
  select(variable, estimate, conf.low, conf.high, p.value, p.value.group, fdr.group) %>% 
  left_join(annotation$nightingale, by = c("variable" = "name.short")) %>% 
  filter(type == "Particle concentration") %>% 
  
  # Reverse the order (add a minus sign)
  mutate(variable = fct_reorder(variable, -name.order)) %>% 
  
  # Before this point: dplyr specifics (data wrangling)
  
  # Below this point: ggplot2 specifics 
  # Note also that I use "+" and not "%>%" after this point
  
  ggplot(aes(x = variable, y = estimate)) + 
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) + 
  coord_flip()

Much better. But did you notice the redundancy in the variable names? They all contain an unnecessary -P. Let’s use the names.pretty instead.

results_scaled_group %>% 
  filter(adj == "Age, gender and weight change") %>% 
  select(variable, estimate, conf.low, conf.high, p.value, p.value.group, fdr.group) %>% 
  left_join(annotation$nightingale, by = c("variable" = "name.short")) %>% 
  filter(type == "Particle concentration") %>% 
  
  # Do the same, but use the names.pretty variable instead
  mutate(name.pretty = fct_reorder(name.pretty, -name.order)) %>% 
  
  # Before this point: dplyr specifics (data wrangling)
  
  # Below this point: ggplot2 specifics 
  # Note also that I use "+" and not "%>%" after this point
  
  ggplot(aes(x = name.pretty, y = estimate)) + 
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) + 
  coord_flip()

Next, we might want to highlight the significance level, here called p.value.group. Let’s see which unique values this variable has.

Colors (default)

unique(results_scaled_group$p.value.group)

## [1] = 0.05  < 0.01  < 0.05  < 0.001
## Levels: = 0.05 < 0.05 < 0.01 < 0.001

This is great. That’s exactly what we want to color by. Add that to the aeshetics to the ggplot call.

results_scaled_group %>% 
  filter(adj == "Age, gender and weight change") %>% 
  select(variable, estimate, conf.low, conf.high, p.value, p.value.group, fdr.group) %>% 
  left_join(annotation$nightingale, by = c("variable" = "name.short")) %>% 
  filter(type == "Particle concentration") %>% 
  mutate(name.pretty = fct_reorder(name.pretty, -name.order)) %>% 
  
  ggplot(aes(x = name.pretty, y = estimate, 
             
             # Color by P value group (significance level)
             color = p.value.group)) + 
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) + 
  coord_flip()

Very cool.

Colors (modified #1)

But what if we are unsatisfied by the default colors that ggplot2 supplied? Then we can choose other colors. I usually pick colors from the palettes offered by the RColorBrewer package. To display them, run this line.

RColorBrewer::display.brewer.all(colorblindFriendly = TRUE)

Here we can color many types of continuous (sequential or diverging) or discrete distributions. We have four groups, and will therefore choose a discrete distribution. I like “Dark2”.

So how do we use this knowledge? I usually modify aesthetics using the "scale_*" family of function. For our particular plot, we want to modify the color aesthetic, and will use the function called "scale_color_brewer)

results_scaled_group %>% 
  filter(adj == "Age, gender and weight change") %>% 
  select(variable, estimate, conf.low, conf.high, p.value, p.value.group, fdr.group) %>% 
  left_join(annotation$nightingale, by = c("variable" = "name.short")) %>% 
  filter(type == "Particle concentration") %>% 
  mutate(name.pretty = fct_reorder(name.pretty, -name.order)) %>% 
  
  ggplot(aes(x = name.pretty, y = estimate, 
             color = p.value.group)) + 
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) + 
  
  # Modify the color aesthetic
  scale_color_brewer(name = "P value", palette = "Dark2") + 
  coord_flip()

We could also choose exactly which colors we want. To do this, we need to retrieve the actual hex codes that make up the palettes.

Say for example that we want a palette that reflects a traffic light for ‘increasing significance’ (red - yellow - green). In addition, we could use a neutral color for non-significance. Let’s find those colors.

RColorBrewer::display.brewer.all()

For “Set1”, there seems to be something that resembles a traffic light. Let’s get the hex codes for those colors.

RColorBrewer::brewer.pal(n = 9, name = "Set1")

## [1] "#E41A1C" "#377EB8" "#4DAF4A" "#984EA3" "#FF7F00" "#FFFF33" "#A65628"
## [8] "#F781BF" "#999999"

Right - that’s what the hex codes look like. Note that the grey color at the end looks a bit different than the rest. So to select our colors of choice, we can simply subset them. To subset, we can use square-brackets, and supply it a vector of the numerical indices that we wish to retain. We would like to retain:

• Grey (index 9)
• Red (index 1)
• Yellow (index 6)
• Green (index 3)

So our code chunk will then look like this:

RColorBrewer::brewer.pal(n = 9, name = "Set1")[c(9, 1, 6, 3)]

## [1] "#999999" "#E41A1C" "#FFFF33" "#4DAF4A"

Now we have our colors. But how to map those onto the correct P value level? In two steps:

1. Supply the P value levels as names for the hex codes, and
2. Use that vector as values in scale_color_manual (see below)

RColorBrewer::brewer.pal(n = 9, name = "Set1")[c(9, 1, 6, 3)] %>% 
  set_names(c(levels(results_scaled_group$p.value.group)))

##    = 0.05    < 0.05    < 0.01   < 0.001 
## "#999999" "#E41A1C" "#FFFF33" "#4DAF4A"

Great - let’s save it to a useful object.

mycols_traffic <- RColorBrewer::brewer.pal(n = 9, name = "Set1")[c(9, 1, 6, 3)] %>% 
  set_names(c(levels(results_scaled_group$p.value.group)))

mycols_traffic

##    = 0.05    < 0.05    < 0.01   < 0.001 
## "#999999" "#E41A1C" "#FFFF33" "#4DAF4A"

And then supply that vector to scale_color_manual.

Colors (modified #2)

results_scaled_group %>% 
  filter(adj == "Age, gender and weight change") %>% 
  select(variable, estimate, conf.low, conf.high, p.value, p.value.group, fdr.group) %>% 
  left_join(annotation$nightingale, by = c("variable" = "name.short")) %>% 
  filter(type == "Particle concentration") %>% 
  mutate(name.pretty = fct_reorder(name.pretty, -name.order)) %>% 
  
  ggplot(aes(x = name.pretty, y = estimate, 
             color = p.value.group)) + 
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) + 
  
  # Switch to scale_color_manual, and set the correct values
  scale_color_manual(name = "P value", values = mycols_traffic) + 
  coord_flip()

In fact, I want to color the inside of the points only, and leave the confidence intervals black. Then I need to do a few more tweaks:

1. Specify the shape of the geom_pointrange to a shape that has both a border and fill, such as shape = 21.
2. Change all color aesthetics to fill aesthetics.

Change from color to fill

results_scaled_group %>% 
  filter(adj == "Age, gender and weight change") %>% 
  select(variable, estimate, conf.low, conf.high, p.value, p.value.group, fdr.group) %>% 
  left_join(annotation$nightingale, by = c("variable" = "name.short")) %>% 
  filter(type == "Particle concentration") %>% 
  mutate(name.pretty = fct_reorder(name.pretty, -name.order)) %>% 
  
  ggplot(aes(x = name.pretty, y = estimate, 
             
             # Change "color" to "fill"
             fill = p.value.group)) + 
  
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), 
                  
                  # Change to a point shape that has both a border and a filling
                  shape = 21) + 
  
  # Switch to scale_fill_manual
  scale_fill_manual(name = "P value", values = mycols_traffic) + 
  coord_flip()

Now most of the actual data-related modifications are in place. However, there are certain other design-related features that we might want to adjust:

1. We can set the theme to a more publication-ready theme
2. We can adjust the x and y labels (Do we need a y label?), and add titles and captions

Theme and labels

results_scaled_group %>% 
  filter(adj == "Age, gender and weight change") %>% 
  select(variable, estimate, conf.low, conf.high, p.value, p.value.group, fdr.group) %>% 
  left_join(annotation$nightingale, by = c("variable" = "name.short")) %>% 
  filter(type == "Particle concentration") %>% 
  mutate(name.pretty = fct_reorder(name.pretty, -name.order)) %>% 
  
  ggplot(aes(x = name.pretty, y = estimate, fill = p.value.group)) + 
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), shape = 21) + 
  scale_fill_manual(name = "P value", values = mycols_traffic) + 
  coord_flip() + 
  
  # Use the classic theme, which is a bit more publication-ready
  theme_classic() + 
  
  # Rename the y-axis label and remove the x-axis label (it doesn't really provide anything)
  labs(y = "Estimate (95 % CI)", x = NULL, 
       
       # Add a title
       title = "Particle concentration")

To continue the polishing efforts, let’s do a few more things:

1. Add a horizontal reference line at estimate = 0 (remember, since we flipped the plotting axes, horizontal lines are vertical)
2. Put the legend at the bottom of the plot
3. Adjust the plot size

results_scaled_group %>% 
  filter(adj == "Age, gender and weight change") %>% 
  select(variable, estimate, conf.low, conf.high, p.value, p.value.group, fdr.group) %>% 
  left_join(annotation$nightingale, by = c("variable" = "name.short")) %>% 
  filter(type == "Particle concentration") %>% 
  mutate(name.pretty = fct_reorder(name.pretty, -name.order)) %>% 
  
  ggplot(aes(x = name.pretty, y = estimate, fill = p.value.group)) + 
  
  # Add a horizontal reference line
  geom_hline(yintercept = 0, color = "grey60", linetype = "dashed") + 
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), shape = 21) + 
  scale_fill_manual(name = "P value", values = mycols_traffic) + 
  coord_flip() + 
  theme_classic() + 
  
  # Move the legend to the bottom
  theme(legend.position = "bottom") + 
  labs(y = "Estimate (95 % CI)", x = NULL, title = "Particle concentration")

When the plot area is this narrow, the legend stretches on the outside of the plot. One solution to this issue is to put the legend icons on two rows.

Legend on two rows

results_scaled_group %>% 
  filter(adj == "Age, gender and weight change") %>% 
  select(variable, estimate, conf.low, conf.high, p.value, p.value.group, fdr.group) %>% 
  left_join(annotation$nightingale, by = c("variable" = "name.short")) %>% 
  filter(type == "Particle concentration") %>% 
  mutate(name.pretty = fct_reorder(name.pretty, -name.order)) %>% 
  
  ggplot(aes(x = name.pretty, y = estimate, fill = p.value.group)) + 
  geom_hline(yintercept = 0, color = "grey60", linetype = "dashed") + 
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), shape = 21) + 
  scale_fill_manual(name = "P value", values = mycols_traffic) + 
  coord_flip() + 
  theme_classic() + 
  theme(legend.position = "bottom") + 
  
  # Put the legend (for the fill aesthetic) on two rows
  guides(fill = guide_legend(nrow = 2)) + 
  labs(y = "Estimate (95 % CI)", x = NULL, title = "Particle concentration")

Finally! I would now say that this plot is complete.

results_scaled_group %>% 
  filter(adj == "Age, gender and weight change") %>% 
  select(variable, estimate, conf.low, conf.high, p.value, p.value.group, fdr.group) %>% 
  left_join(annotation$nightingale, by = c("variable" = "name.short")) %>% 
  
  # Filter on more variable types (defined in the annotation file)
  filter(type %in% c("Particle concentration", "Apolipoproteins") |
           lipid.subclasses == "no" & type %in% c(
             "Cholesterol", "Triglycerides", "Phospholipids")
         )

First try

results_scaled_group %>% 
  filter(adj == "Age, gender and weight change") %>% 
  select(variable, estimate, conf.low, conf.high, p.value, p.value.group, fdr.group) %>% 
  left_join(annotation$nightingale, by = c("variable" = "name.short")) %>% 
  
  # Filter on more variable types (defined in the annotation file)
  filter(type %in% c("Particle concentration", "Apolipoproteins") |
           lipid.subclasses == "no" & type %in% c("Cholesterol", 
                                                  "Triglycerides", 
                                                  "Phospholipids")) %>% 
  mutate(name.pretty = fct_reorder(name.pretty, -name.order)) %>% 
  
  ggplot(aes(x = name.pretty, y = estimate, fill = p.value.group)) + 
  geom_hline(yintercept = 0, color = "grey60", linetype = "dashed") + 
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), shape = 21) + 
  scale_fill_manual(name = "P value", values = mycols_traffic) + 
  coord_flip() + 
  theme_classic() + 
  theme(legend.position = "bottom") + 
  guides(fill = guide_legend(nrow = 2)) + 
  labs(y = "Estimate (95 % CI)", x = NULL)

Alright. The figure is now higher (note that I doubled the fig.height in the code chunk header!).

But something is wrong. We need to group the variables based on the types. To this, we use facets.

Facet by variable types

results_scaled_group %>% 
  filter(adj == "Age, gender and weight change") %>% 
  select(variable, estimate, conf.low, conf.high, p.value, p.value.group, fdr.group) %>% 
  left_join(annotation$nightingale, by = c("variable" = "name.short")) %>% 
  filter(type %in% c("Particle concentration", "Apolipoproteins") |
           lipid.subclasses == "no" & type %in% c("Cholesterol", 
                                                  "Triglycerides", 
                                                  "Phospholipids")) %>% 
  mutate(name.pretty = fct_reorder(name.pretty, -name.order)) %>% 
  
  ggplot(aes(x = name.pretty, y = estimate, fill = p.value.group)) + 
  geom_hline(yintercept = 0, color = "grey60", linetype = "dashed") + 
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), shape = 21) + 
  scale_fill_manual(name = "P value", values = mycols_traffic) + 
  coord_flip() + 
  
  # Create a facet grid by variable types in the rows
  facet_grid(rows = vars(type)) + 
  theme_classic() + 
  theme(legend.position = "bottom") + 
  guides(fill = guide_legend(nrow = 2)) + 
  labs(y = "Estimate (95 % CI)", x = NULL)

This is obviously not what we want. We need to free up the scales and plotting space within each facet.

results_scaled_group %>% 
  filter(adj == "Age, gender and weight change") %>% 
  select(variable, estimate, conf.low, conf.high, p.value, p.value.group, fdr.group) %>% 
  left_join(annotation$nightingale, by = c("variable" = "name.short")) %>% 
  filter(type %in% c("Particle concentration", "Apolipoproteins") |
           lipid.subclasses == "no" & type %in% c("Cholesterol", 
                                                  "Triglycerides", 
                                                  "Phospholipids")) %>% 
  mutate(name.pretty = fct_reorder(name.pretty, -name.order)) %>% 
  
  ggplot(aes(x = name.pretty, y = estimate, fill = p.value.group)) + 
  geom_hline(yintercept = 0, color = "grey60", linetype = "dashed") + 
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), shape = 21) + 
  scale_fill_manual(name = "P value", values = mycols_traffic) + 
  coord_flip() + 
  facet_grid(rows = vars(type), 
             
             # Free up the scales and plotting space within each facet
             scales = "free", space = "free") + 
  theme_classic() + 
  theme(legend.position = "bottom") + 
  guides(fill = guide_legend(nrow = 2)) + 
  labs(y = "Estimate (95 % CI)", x = NULL)

I don’t completely like the ugly black border around the strips. Also, I want to be able to read the strip text without tilting my head (let’s rotate them!). Again we modify non-data ink using themes.

results_scaled_group %>% 
  filter(adj == "Age, gender and weight change") %>% 
  select(variable, estimate, conf.low, conf.high, p.value, p.value.group, fdr.group) %>% 
  left_join(annotation$nightingale, by = c("variable" = "name.short")) %>% 
  filter(type %in% c("Particle concentration", "Apolipoproteins") |
           lipid.subclasses == "no" & type %in% c("Cholesterol", 
                                                  "Triglycerides", 
                                                  "Phospholipids")) %>% 
  mutate(name.pretty = fct_reorder(name.pretty, -name.order)) %>% 
  
  ggplot(aes(x = name.pretty, y = estimate, fill = p.value.group)) + 
  geom_hline(yintercept = 0, color = "grey60", linetype = "dashed") + 
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), shape = 21) + 
  scale_fill_manual(name = "P value", values = mycols_traffic) + 
  coord_flip() + 
  facet_grid(rows = vars(type), scales = "free", space = "free") + 
  theme_classic() + 
  theme(legend.position = "bottom", 
        
        # Remove ugly black border around strips
        strip.background.y = element_blank(), 
        
        # Rotate the strip text so that it's easier to read
        strip.text.y = element_text(angle = 0)) + 
  guides(fill = guide_legend(nrow = 2)) + 
  labs(y = "Estimate (95 % CI)", x = NULL)

This, unfortunately, makes the figure a bit too wide. There is a lot of whitespace on the right hand side. There are several ways to deal with this, one of which is to simply shorten the faceting names. I’ve already prepared such a variable, so let’s facet by that one instead.

results_scaled_group %>% 
  filter(adj == "Age, gender and weight change") %>% 
  select(variable, estimate, conf.low, conf.high, p.value, p.value.group, fdr.group) %>% 
  left_join(annotation$nightingale, by = c("variable" = "name.short")) %>% 
  filter(type %in% c("Particle concentration", "Apolipoproteins") |
           lipid.subclasses == "no" & type %in% c("Cholesterol", 
                                                  "Triglycerides", 
                                                  "Phospholipids")) %>% 
  mutate(name.pretty = fct_reorder(name.pretty, -name.order)) %>% 
  
  ggplot(aes(x = name.pretty, y = estimate, fill = p.value.group)) + 
  geom_hline(yintercept = 0, color = "grey60", linetype = "dashed") + 
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), shape = 21) + 
  scale_fill_manual(name = "P value", values = mycols_traffic) + 
  coord_flip() + 
  
  # Change to "type.short", which displays nicer on the facets
  facet_grid(rows = vars(type.short), scales = "free", space = "free") + 
  theme_classic() + 
  theme(legend.position = "bottom", 
        strip.background.y = element_blank(), 
        strip.text.y = element_text(angle = 0)) + 
  guides(fill = guide_legend(nrow = 2)) + 
  labs(y = "Estimate (95 % CI)", x = NULL)

That’s better. We can just add the whole facet/group name in an abbreviation list.

In a final tweak, I want to manually switch the order in which a few of the variables are presented. For example, I often want to show ApoB before ApoA1, and I want to move the variables EstC, FreeC, and Remnant-C above the specific particles, within the Cs/Cholesterol facet. We can do all this before we plot the figure, inside the call to mutate.

results_scaled_group %>% 
  filter(adj == "Age, gender and weight change") %>% 
  select(variable, estimate, conf.low, conf.high, p.value, p.value.group, fdr.group) %>% 
  left_join(annotation$nightingale, by = c("variable" = "name.short")) %>% 
  filter(type %in% c("Particle concentration", "Apolipoproteins") |
           lipid.subclasses == "no" & type %in% c("Cholesterol", 
                                                  "Triglycerides", 
                                                  "Phospholipids")) %>% 
  mutate(name.pretty = fct_reorder(name.pretty, -name.order) %>% 
           
           # Relevel a few specific variables
           fct_relevel( 
             
             # Pull ApoB to the start of it's panel
             "ApoB", 
             
             # Pull *EstC*, *FreeC*, and *Remnant-C* to the start of their panel (after "Serum-C")
             "Remnant-C", "FreeC", "EstC", "Serum-C", 
             
             after = Inf
             
             )) %>% 
  
  ggplot(aes(x = name.pretty, y = estimate, fill = p.value.group)) + 
  geom_hline(yintercept = 0, color = "grey60", linetype = "dashed") + 
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), shape = 21) + 
  scale_fill_manual(name = "P value", values = mycols_traffic) + 
  coord_flip() + 
  facet_grid(rows = vars(type.short), scales = "free", space = "free") + 
  theme_classic() + 
  theme(legend.position = "bottom", 
        strip.background.y = element_blank(), 
        strip.text.y = element_text(angle = 0)) + 
  guides(fill = guide_legend(nrow = 2)) + 
  labs(y = "Estimate (95 % CI)", x = NULL)

And THAT … I consider a nice plot. Here is the final code snippet.

results_scaled_group %>% 
  filter(adj == "Age, gender and weight change") %>% 
  select(variable, estimate, conf.low, conf.high, p.value, p.value.group, fdr.group) %>% 
  left_join(annotation$nightingale, by = c("variable" = "name.short")) %>% 
  filter(type %in% c("Particle concentration", "Apolipoproteins") |
           lipid.subclasses == "no" & type %in% c("Cholesterol", "Triglycerides", "Phospholipids")) %>% 
  mutate(name.pretty = fct_reorder(name.pretty, -name.order) %>% 
           fct_relevel("ApoB", "Remnant-C", "FreeC", "EstC", "Serum-C", after = Inf)) %>% 
  ggplot(aes(x = name.pretty, y = estimate, fill = p.value.group)) + 
  geom_hline(yintercept = 0, color = "grey60", linetype = "dashed") + 
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), shape = 21) + 
  scale_fill_manual(name = "P value", values = mycols_traffic) + 
  coord_flip() + 
  facet_grid(rows = vars(type.short), scales = "free", space = "free") + 
  theme_classic() + 
  theme(legend.position = "bottom", 
        strip.background.y = element_blank(), 
        strip.text.y = element_text(angle = 0)) + 
  guides(fill = guide_legend(nrow = 2)) + 
  labs(y = "Estimate (95 % CI)", x = NULL)

Save the final figure

Let’s save this figure using ggsave. This function remembers the last plot that you produced (here, in the RMarkdown document), and then saves that to a file that you specify.

Use this code to save PNG (raster graphics) files.

ggsave(filename = "../results/figures/forest-lipid-classes.png", 
       width = 5, height = 8, dpi = 300)

Use this code to save PDF (vector graphics) files.

ggsave(filename = "../results/figures/forest-lipid-classes.pdf", 
       width = 5, height = 8, device = cairo_pdf)

Note that I set the device to cairo_pdf from the Cairo package. This is useful when we save PDF files, because it can render special characters, like the more-than-or-equal sign (>=, ≥) that I use to signify more-than-or-equal-to 0.05 (for P values) or 0.20 (for FDR values). If I use the regular pdf device, this would look like a equal sign (=), which is obviously incorrect.